Exercise 1: Create study area map using “ggplot”

All datasets are found in the “data”-folder.

Here, we use the package ‘rnaturalearth’ to download countries in a ESRI shapefile format. These are available for download here: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-0-countries/

Coordinate reference systems (CRS)
CRS provide a standardized way of describing locations to describe geographic data. The CRS that is chosen depends on when the data was collected, the purpose of the data, etc. It is necessary to transform vector and raster data to a common CRS when data with different CRS are combined.

Setting CRS to SWEREF99 (Sweden). CRS can be referenced by its:

  1. EPSG code, CRS(“+init=epsg:3006”) (see http://www.epsg-registry.org/ and http://spatialreference.org/) or by

  2. proj4spring; “+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs”

There are two general options:

  1. unprojected (a.k.a. Geographic): Latitude/Longitude for referencing location on the ellipsoid Earth (wgs84), and

  2. projected: Easting/Northing for referencing location on 2D representations of Earth (the creation of maps) e.g., SWEREF99.

More reading:
https://www.nceas.ucsb.edu/sites/default/files/2020-04/OverviewCoordinateReferenceSystems.pdf

Load vector data

How to download using “rnaturalearth” and save the shp locally

I will open the saved .shp for the exercise:

## load country boarders -----------


#' downladed and saved from rnaturalearth (above)
swe_no_sweref <-  input_vector %>% 

  # value = T: return vector containing the matching elements
  grep(pattern = "swe_no", value = T) %>% 

  # read ESRI Shapefile object 
  sf::st_read() %>% 

  #convert geometry object into an sf object
  st_as_sf() 
## Reading layer `swe_no' from data source 
##   `/Users/hirn0001/Library/CloudStorage/OneDrive-Sverigeslantbruksuniversitet/Undervisning/Teaching and courses/GS_VMAS Data handling and illustrations/2025/create-maps/input_data/vector/country/swe_no.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -362704.5 ymin: -6097519 xmax: 1114145 ymax: 8972413
## Projected CRS: SWEREF99 TM
# check coordinate system 
st_crs(swe_no_sweref)$proj4string
## [1] "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
# Create a bounding box as an sf object for layers ----
b_box <- c(xmin = 350432, ymin = 7179284, 
           xmax = 857780, ymax =7444230)

# Now intersect the full world's coastline with the bounding box 
suppressWarnings({ 
sweden.c <- st_crop(swe_no_sweref, b_box)
plot(sweden.c)
}) 

class(sweden.c)
## [1] "sf"         "data.frame"
### water data ----

# import lakes
lakes <-  input_vector %>%
  grep(pattern = "lakes", value = T) %>%
  sf::st_read() %>%
  st_as_sf() 
## Reading layer `ne_10m_lakes' from data source 
##   `/Users/hirn0001/Library/CloudStorage/OneDrive-Sverigeslantbruksuniversitet/Undervisning/Teaching and courses/GS_VMAS Data handling and illustrations/2025/create-maps/input_data/vector/ne_10m_lakes/ne_10m_lakes.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1355 features and 41 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -165.9656 ymin: -50.66967 xmax: 177.1544 ymax: 81.95521
## Geodetic CRS:  WGS 84
class(lakes)
## [1] "sf"         "data.frame"
st_crs(lakes)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
# Transform or convert coordinates of simple feature (sf)
lakes <- st_transform(lakes, crs = 3006) 
suppressWarnings({ 
  
# crop an sf object to a specific rectangle (here based on extent of "sweden.c"-shp)
# sf_use_s2(FALSE) # use if errors occur below
lakes <- st_crop(lakes, sweden.c) 

class(lakes)
plot(lakes)
}) 

Add points of interest (cities)

Create xy points with labels. Here creating sf-object and setting the crs.

# Create points of interest (coordinates from QGIS or google maps for demonstration)
places <- data.frame(ID= c("Arjeplog", "Arvidsjaur", "Sorsele"),
                     y = c(7328407.945, 7280856.153, 7270326.684),
                     x = c(630040.101, 693262.414, 617477.638)) %>%
  st_as_sf(coords = c("x", "y"),  # create sf object 
           crs=st_crs(sweden.c))  # set same crs as for sweden.c     
class(places)
## [1] "sf"         "data.frame"
# st_crs(places) # check projection 

Make plot using ggplot

library(ggspatial)

ggplot() +
  
    # plot sweden and lakes 
    geom_sf(aes(geometry = sweden.c$geometry)) + # note: geom_sf() for vector data 
  
    geom_sf(data = sweden.c,                     # plot sweden border vector
          aes(geometry = sweden.c$geometry),
          fill= "antiquewhite")+

      geom_sf(data = lakes,                      # plot lakes vector
              fill="lightblue") +

  # layer_spatial(data=elevation) +
  
  # add cities 
  geom_sf(data = places) +
  
  # add city names 
  geom_sf_text(data = places, aes(label = ID),
               vjust=-0.5)+
  
  # Add text on map
  annotate(geom = "text", 
             x = 484609, y = 7404230, 
             label = "Norway", 
             fontface = "italic", 
             color = "darkgray", 
           size = 6) +
  
    annotate(geom = "text", 
             x = 784609, y = 7404230, 
             label = "Sweden", 
             fontface = "italic", 
             color = "darkgray", 
           size = 6) +
  
  # Add north arrow
  annotation_north_arrow(location = "bl", 
                       which_north = "true", 
                         height = unit(1, "cm"),
  width = unit(1, "cm"),
                       pad_x = unit(0.5, "cm"), # horizontal align
                       pad_y = unit(8, "cm"), # vertical align
                       style = north_arrow_fancy_orienteering) +
  #add annotation scale
  annotation_scale(location = "bl", 
                     width_hint = 0.5) +
  
  
  # keep new projection for ggplot
  coord_sf(crs = st_crs(3006), 
  datum = sf::st_crs(3006),             
  expand = FALSE) +
  
  # add grid lines 
  theme(panel.grid.major = element_line(color = gray(.5), 
                                        linetype = "dashed", 
                                        linewidth  = 0.1), 
        panel.background = element_rect(fill = "lightblue")) 

### examplefor removing or adding titles 
# theme(axis.title.x = element_blank(),    # remove 
#       axis.title.y = element_blank()) + 
  # xlab("Longitude") + ylab("Latitude")   # add

ggsave("output_data/study_area.png", 
       width = 8, 
       height = 6, 
       dpi = 400)

Optional: Create a polygon and add to map

Exercise 2 interactive maps

reindeer_data <- readRDS("input_data/vector/animal_data/reindeer_data.rds") %>%
   dplyr::select(Collar_ID,  id, year, x_, y_, t_, group) 

str(reindeer_data)
## 'data.frame':    5951 obs. of  7 variables:
##  $ Collar_ID: chr  "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
##  $ id       : chr  "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
##  $ year     : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ x_       : num  733706 733653 733836 733828 733683 ...
##  $ y_       : num  7314782 7314737 7314807 7314851 7314738 ...
##  $ t_       : chr  "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
##  $ group    : chr  "Fed" "Fed" "Fed" "Fed" ...
# character to factor needed for grouping colors
reindeer_data$group <- as.factor(reindeer_data$group) 
unique(reindeer_data$Collar_ID)
##  [1] "FIT_421952" "FIT_421943" "FIT_421941" "FIT_421940" "FIT_421949"
##  [6] "FIT_421910" "FIT_421939" "FIT_421946" "FIT_421906" "FIT_421956"
## [11] "FIT_421944" "FIT_421937" "FIT_421948" "FIT_421908" "FIT_421935"
## [16] "FIT_421942" "FIT_421954" "FIT_421907" "FIT_421932" "FIT_421931"
## [21] "FIT_421934" "FIT_421945" "FIT_421938" "FIT_421947" "FIT_421936"
## [26] "FIT_421953" "FIT_421909" "FIT_421950" "FIT_421951" "FIT_421933"
library(mapview)
mapid <- reindeer_data %>%
   filter(Collar_ID =="FIT_421952" | Collar_ID =="FIT_421943"| Collar_ID =="FIT_421941")  %>%
  st_as_sf(coords = c("x_", "y_"), crs=3006)

AOI <- grep("modified_new.shp", input_vector, value = T) %>%
  sf::read_sf() %>%
  sf::st_transform(crs = 3006)

pal <-  mapviewPalette("mapviewVectorColors") # optional: mapviewTopoColors, mapviewRasterColors, mapviewVectorColors, mapviewSpectralColors
mapid %>%
  mapview(zcol="Collar_ID", 
          cex = "Collar_ID", 
          col.regions = pal(3)) + 
  mapview(
    AOI,
    alpha.regions = 0,   # no fill
    color = "black",     # border color
    lwd = 2              # border width
  ) # or simply add + AOI 
### Some more examples....



mapid2 <- reindeer_data %>%
  filter(Collar_ID %in% c("FIT_421952","FIT_421943","FIT_421941")) %>%
  mutate(t_ = as.POSIXct(t_, format = "%Y-%m-%d %H:%M:%S", tz = "UTC")) %>%
  arrange(Collar_ID, t_) %>%
  st_as_sf(coords = c("x_", "y_"), crs = 3006)

tracks <- mapid2 %>%
  group_by(Collar_ID) %>%
 summarise(do_union = FALSE) %>%     # keeps point order within groups
  st_cast("LINESTRING")


mapview(tracks, 
        zcol = "Collar_ID", 
        lwd = 3, 
        layer.name = "Tracks") +
  mapview(mapid, zcol = "Collar_ID", cex = 3, layer.name = "Fixes") 

Exercise 3 (optional): Maps in R

Pre-define crs using ESPG-code or proj4string by:

# projection to use as proj4string 
crs_to_use <- "+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs" 

Load background data

Now let’s start by loading the maps using the “raster”-package. Another option: “terra”-package using rast()-function)

## [1] "input_data/raster/PV_2021_cc.tif"    
## [2] "input_data/raster/PV_2021_dem.tif"   
## [3] "input_data/raster/PV_2021_houses.tif"
## [4] "input_data/raster/PV_2021_lichen.tif"
## [5] "input_data/raster/PV_2021_nmd.tif"   
## [6] "input_data/raster/PV_2021_roads.tif" 
## [7] "input_data/raster/PV_2021_slope.tif" 
## [8] "input_data/raster/PV_2021_TRI.tif"
## [1] "Clear_cuts_2021"            "PV_2021_dem"               
## [3] "prox10_houses_stakke"       "lichenmap_stakke_original" 
## [5] "nmd10_stakke_ungeneralized" "prox10_stakke_resample"    
## [7] "slope10m_QGIS"              "ruggedness10_QGIS"
## [1] "+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs"

## Reading layer `spatial_files_availability_PV_winter_2021_shorter_modified_new' from data source `/Users/hirn0001/Library/CloudStorage/OneDrive-Sverigeslantbruksuniversitet/Undervisning/Teaching and courses/GS_VMAS Data handling and illustrations/2025/create-maps/input_data/vector/HR/spatial_files_availability_PV_winter_2021_shorter_modified_new.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 716285.4 ymin: 7313268 xmax: 739293.8 ymax: 7321582
## Projected CRS: SWEREF99 TM
## [1] "+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs"

Add animal data (GPS points)

Load prepared reindeer data

## 'data.frame':    5951 obs. of  7 variables:
##  $ Collar_ID: chr  "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
##  $ id       : chr  "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
##  $ year     : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ x_       : num  733706 733653 733836 733828 733683 ...
##  $ y_       : num  7314782 7314737 7314807 7314851 7314738 ...
##  $ t_       : chr  "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
##  $ group    : chr  "Fed" "Fed" "Fed" "Fed" ...
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"

Explore and select colors

## [1] "#00A600" "#01A600" "#03A700" "#04A700" "#05A800" "#07A800"

Plot and save all together

## Warning in plot.sf(HR_area, col = NA, border = 1, add = T): ignoring all but
## the first attribute